packages <- c("CIMseq", "CIMseq.data", "tidyverse", "scSeqTools")
purrr::walk(packages, library, character.only = TRUE)
rm(packages)
##DATA
load('../data/CIMseqData.rda')
load('../data/sObj.rda')
#rename classes
renameClasses <- function(class) {
case_when(
class == "0" ~ "C.Stem.proximal",
class == "1" ~ "C.Stem.distal",
class == "2" ~ "SI.Stem",
class == "3" ~ "C.Goblet.distal",
class == "4" ~ "SI.Enterocytes",
class == "5" ~ "SI.Goblet",
class == "6" ~ "C.Goblet.proximal",
class == "7" ~ "Enteroendocrine",
class == "8" ~ "Tufft",
class == "9" ~ "SI.Paneth",
class == "10" ~ "Blood",
TRUE ~ "error"
)
}
cOrder <- c(
"C.Stem.proximal", "C.Goblet.proximal", "C.Stem.distal", "C.Goblet.distal",
"SI.Goblet", "SI.Paneth", "SI.Stem", "SI.Enterocytes",
"Enteroendocrine", "Tufft", "Blood"
)
getData(cObjSng, "classification") <- renameClasses(getData(cObjSng, "classification"))
features <- getData(cObjMul, "features")
names(features) <- renameClasses(names(features))
getData(cObjMul, "features") <- features
fractions <- getData(sObj, "fractions")
colnames(fractions) <- renameClasses(colnames(fractions))
sObj@fractions <- fractions
plotUnsupervisedClass(cObjSng, cObjMul)
plotUnsupervisedMarkers(
cObjSng, cObjMul,
c("Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1"),
pal = RColorBrewer::brewer.pal(8, "Set1")
)
plotUnsupervisedMarkers(
cObjSng, cObjMul, "Hoxb13",
pal = RColorBrewer::brewer.pal(8, "Set1")
)
plotUnsupervisedMarkers(
cObjSng, cObjMul, "Mki67",
pal = RColorBrewer::brewer.pal(8, "Set1")
)
classHeatmap(
data = data.frame(
gene = rownames(getData(cObjMul, "counts"))[getData(cObjMul, "features")],
class = names(getData(cObjMul, "features"))
),
counts.log = getData(cObjSng, "counts.log"),
classes = tibble(
sample = colnames(getData(cObjSng, "counts")),
class = getData(cObjSng, "classification")
)
)
violinPlot <- function(cpm, genes, classes, class) {
cpm[genes, ] %>%
t() %>%
matrix_to_tibble("sample") %>%
gather(gene, cpm, -sample) %>%
full_join(tibble(sample = colnames(cpm), class = classes), by = "sample") %>%
ggplot() +
geom_jitter(aes(class, cpm), height = 0, size = 0.1, alpha = 0.25) +
facet_wrap(~gene, scales = "free_y") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90)) +
labs(title = paste0(class, " genes"))
}
getGenes <- function(
cpm, classes, features, class, ntop = 10
) {
fc <- foldChangePerClass(
cpm[features[names(features) == class], ],
tibble(sample = colnames(cpm), class = classes)
)
rownames(fc[order(fc[, class], decreasing = TRUE), ])[1:ntop]
}
c <- "C.Stem.distal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
c <- "C.Stem.proximal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
c <- "C.Goblet.distal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
c <- "C.Goblet.proximal"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
c <- "SI.Stem"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
c <- "SI.Goblet"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
c <- "SI.Enterocytes"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
c <- "SI.Paneth"
g <- getGenes(getData(cObjSng, "counts.cpm"), getData(cObjSng, "classification"), getData(cObjMul, "features"), c)
violinPlot(getData(cObjSng, "counts.cpm"), g, getData(cObjSng, "classification"), c)
plotSwarmCircos(sObj, cObjSng, cObjMul, weightCut = 0, classOrder = cOrder)
Only detected duplicates and triplicates.
Only ERCC estimated cell number <= 4.
Weight cutoff = 5.
adj <- adjustFractions(cObjSng, cObjMul, sObj, binary = TRUE)
samples <- rownames(adj)
rs <- rowSums(adj)
keep <- rs == 2 | rs == 3
plotSwarmCircos(
filterSwarm(sObj, keep), cObjSng, cObjMul, weightCut = 5,
classOrder = cOrder, theoretical.max = 4
)